Install and/or load packages

# install.packages("tidyverse")
# install.packages("haven")
# install.packages("lme4")
# install.packages("lmerTest")
# install.packages("effects")

library(tidyverse)
library(haven)
library(lme4)
library(lmerTest)
library(effects)

# set default contrasts
options(contrasts = c("contr.sum", "contr.poly"))

Custom functions

lmer_vpc <- function(fit) {
  fit %>% 
    VarCorr() %>% 
    as_tibble() %>% 
    mutate(vpc = vcov / sum(vcov),
           vpc_pct = paste0(round(vpc * 100, 2), "%"))
}

Data

impact_data <- read_spss("data/cleandata/final_imPACT_long_desire_stress_jessica_nick.sav")
impact_data_long <- read_spss("data/cleandata/15_05_14_imPACT_use_for_big_T.sav")

# replace NaN with NA
impact_data[is.na(impact_data)] <- NA
impact_data_long[is.na(impact_data_long)] <- NA

Join with uncentered BMI

impact_data <- impact_data %>% 
  left_join(
    impact_data_long %>% 
  distinct(ID, BMI_B), by = "ID"
  )

Change and/or create new variables

impact_data <- impact_data %>% 
  mutate(gender_chr = case_when(Gender == 0 ~ "Man", Gender == 1 ~ "Woman"),
         gender_fct = factor(gender_chr),
         relationship_status_chr = case_when(relstatcoded == 1 ~ "Single", relstatcoded == 2 ~ "Casual", relstatcoded == 3 ~ "Committed"),
         relationship_status_fct = factor(relationship_status_chr, levels = c("Single", "Casual", "Committed")),
         nicotine_chr = case_when(nicotine == 0 ~ "No", relstatcoded == 1 ~ "Yes"),
         nicotine_fct = factor(nicotine_chr),
         sex_med_chr = case_when(sexmeds == 0 ~ "No", sexmeds == 1 ~ "Yes"),
         sex_med_fct = factor(sex_med_chr),
         testosterone_med_chr = case_when(TMeds == 0 ~ "No", TMeds == 1 ~ "Yes"),
         testosterone_med_fct = factor(testosterone_med_chr),
         cortisol_med_chr = case_when(CMeds == 0 ~ "No", CMeds == 1 ~ "Yes"),
         cortisol_med_fct = factor(cortisol_med_chr),
         hormonal_contraceptive_chr = case_when(HC == 0 ~ "No", HC == 1 ~ "Yes"),
         hormonal_contraceptive_fct = factor(hormonal_contraceptive_chr),
         illness_chr = case_when(illness == 0 ~ "No", illness == 1 ~ "Cold or minor includes unspecified", illness == 2 ~ "More intense or chronic"),
         illness_fct = factor(illness_chr, levels = c("No", "Cold or minor includes unspecified", "More intense or chronic")),
         relationship_status_fct = factor(relationship_status_chr))

Set contrasts

contrasts(impact_data$relationship_status_fct) <- contr.poly(n = 3)
contrasts(impact_data$illness_fct) <- contr.poly(n = 3)
contrasts(impact_data$hormonal_contraceptive_fct) <- c(-1, 1)

Subset data by exclusion criteria and gender

# men and women
impact_data_nodrugs <- impact_data %>% 
  filter(nicotine_chr == "No" & sex_med_chr == "No" & testosterone_med_chr == "No" & cortisol_med_chr == "No") %>% 
  mutate(time_since_waking_dbl = parse_number(as.character(timesincewake)),
         time_since_waking_c = time_since_waking_dbl - mean(time_since_waking_dbl, na.rm = TRUE),
         bmi_c = BMI_B - mean(BMI_B, na.rm = TRUE),
         testosterone_screen_c = Testosteronescreen - mean(Testosteronescreen, na.rm = TRUE),
         cortisol_screen_c = Cortscreen - mean(Cortscreen, na.rm = TRUE),
         pss_c = PSS - mean(PSS, na.rm = TRUE))

# men
impact_data_nodrugs_men <- impact_data %>% 
  filter(gender_chr == "Man" & nicotine_chr == "No" & sex_med_chr == "No" & testosterone_med_chr == "No" & cortisol_med_chr == "No") %>% 
  mutate(time_since_waking_dbl = parse_number(as.character(timesincewake)),
         time_since_waking_c = time_since_waking_dbl - mean(time_since_waking_dbl, na.rm = TRUE),
         bmi_c = BMI_B - mean(BMI_B, na.rm = TRUE),
         testosterone_screen_c = Testosteronescreen - mean(Testosteronescreen, na.rm = TRUE),
         cortisol_screen_c = Cortscreen - mean(Cortscreen, na.rm = TRUE),
         pss_c = PSS - mean(PSS, na.rm = TRUE))

# women
impact_data_nodrugs_women <- impact_data %>% 
  filter(gender_chr == "Woman" & nicotine_chr == "No" & sex_med_chr == "No" & testosterone_med_chr == "No" & cortisol_med_chr == "No") %>% 
  mutate(time_since_waking_dbl = parse_number(as.character(timesincewake)),
         time_since_waking_c = time_since_waking_dbl - mean(time_since_waking_dbl, na.rm = TRUE),
         bmi_c = BMI_B - mean(BMI_B, na.rm = TRUE),
         testosterone_screen_c = Testosteronescreen - mean(Testosteronescreen, na.rm = TRUE),
         cortisol_screen_c = Cortscreen - mean(Cortscreen, na.rm = TRUE),
         pss_c = PSS - mean(PSS, na.rm = TRUE))

Exploratory Plots

Sum Dyadic Sexual Desire (6-62 points) by gender (boxplots)

impact_data_nodrugs %>% 
  ggplot(mapping = aes(x = ordered(Time_AdjZero), y = SDI_dyad, fill = gender_fct)) +
  geom_boxplot() +
  scale_x_discrete(breaks = seq(0, 11, 1)) +
  scale_y_continuous(breaks = seq(0, 70, 10)) +
  scale_fill_manual(values = c("#0072B2", "#CC79A7")) +
  labs(x = "Session (0-11 session numbers)", y = NULL, title = "Sum Dyadic Sexual Desire (6-62 points) by gender", fill = NULL) +
  theme(legend.position = "top")
## Warning: Removed 25 rows containing non-finite values (stat_boxplot).

Sum Dyadic Sexual Desire (6-62 points) by gender (cubic splines)

impact_data_nodrugs %>% 
  ggplot(mapping = aes(x = Time_AdjZero, y = SDI_dyad, fill = gender_fct, color = gender_fct)) +
  geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs", k = 5), se = FALSE) +
  geom_point(alpha = 0.25, size = 3) +
  scale_x_continuous(breaks = seq(0, 11, 1)) +
  scale_y_continuous(breaks = seq(0, 70, 10)) +
  scale_fill_manual(values = c("#0072B2", "#CC79A7")) +
  scale_color_manual(values = c("#0072B2", "#CC79A7")) +
  labs(x = "Session (0-11 session numbers)", y = NULL, title = "Sum Dyadic Sexual Desire (6-62 points) by gender", subtitle = "Cubic Spline Smooth, k = 5", fill = NULL, color = NULL) +
  theme(legend.position = "top")
## Warning: Removed 25 rows containing non-finite values (stat_smooth).
## Warning: Removed 25 rows containing missing values (geom_point).

Sum Solitary Sexual Desire (2-23 points) by gender (boxplots)

impact_data_nodrugs %>% 
  ggplot(mapping = aes(x = ordered(Time_AdjZero), y = SDI_sol, fill = gender_fct)) +
  geom_boxplot() +
  scale_x_discrete(breaks = seq(0, 11, 1)) +
  scale_y_continuous(breaks = seq(0, 25, 5)) +
  scale_fill_manual(values = c("#0072B2", "#CC79A7")) +
  labs(x = "Session (0-11 session numbers)", y = NULL, title = "Sum Solitary Sexual Desire (2-23 points) by gender", fill = NULL) +
  theme(legend.position = "top")
## Warning: Removed 23 rows containing non-finite values (stat_boxplot).

Sum Solitary Sexual Desire (2-23 points) by gender (cubic splines)

impact_data_nodrugs %>% 
  ggplot(mapping = aes(x = Time_AdjZero, y = SDI_sol, fill = gender_fct, color = gender_fct)) +
  geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs", k = 5), se = FALSE) +
  geom_point(alpha = 0.25, size = 3) +
  scale_x_continuous(breaks = seq(0, 11, 1)) +
  scale_y_continuous(breaks = seq(0, 25, 5)) +
  scale_fill_manual(values = c("#0072B2", "#CC79A7")) +
  scale_color_manual(values = c("#0072B2", "#CC79A7")) +
  labs(x = "Session (0-11 session numbers)", y = NULL, title = "Sum Solitary Sexual Desire (2-23 points) by gender", subtitle = "Cubic Spline Smooth, k = 5", fill = NULL, color = NULL) +
  theme(legend.position = "top")
## Warning: Removed 23 rows containing non-finite values (stat_smooth).
## Warning: Removed 23 rows containing missing values (geom_point).

Relationship among testosterone, cortisol, and dyadic sexual desire in men

impact_data_nodrugs_men %>%
  mutate(cortisol_screen_c_q = cut(cortisol_screen_c, breaks = quantile(cortisol_screen_c, probs = seq(0, 1, 1 / 5), na.rm = TRUE), include.lowest = TRUE)) %>% 
  filter(!is.na(cortisol_screen_c_q)) %>% 
  ggplot(mapping = aes(x = testosterone_screen_c, y = SDI_dyad)) +
  geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"), se = FALSE, color = "#0072B2") +
  geom_point(alpha = 0.25, size = 3, color = "#0072B2") +
  facet_wrap(facets = ~ cortisol_screen_c_q) +
  scale_y_continuous(breaks = seq(0, 70, 10)) +
  labs(x = "Testosterone (mean-centered)", y = NULL, title = "Sum Dyadic Sexual Desire (6-62 points) by mean-centered cortisol quintile in men") +
  theme(legend.position = "top")
## Warning: Removed 17 rows containing non-finite values (stat_smooth).
## Warning: Removed 17 rows containing missing values (geom_point).

Relationship among testosterone, cortisol, and solitary sexual desire in men

impact_data_nodrugs_men %>%
  mutate(cortisol_screen_c_q = cut(cortisol_screen_c, breaks = quantile(cortisol_screen_c, probs = seq(0, 1, 1 / 5), na.rm = TRUE), include.lowest = TRUE)) %>% 
  filter(!is.na(cortisol_screen_c_q)) %>% 
  ggplot(mapping = aes(x = testosterone_screen_c, y = SDI_sol)) +
  geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"), se = FALSE, color = "#0072B2") +
  geom_point(alpha = 0.25, size = 3, color = "#0072B2") +
  facet_wrap(facets = ~ cortisol_screen_c_q) +
  scale_y_continuous(breaks = seq(0, 25, 5)) +
  labs(x = "Testosterone (mean-centered)", y = NULL, title = "Sum Solitary Sexual Desire (2-23 points) by mean-centered cortisol quintile in men") +
  theme(legend.position = "top")
## Warning: Removed 20 rows containing non-finite values (stat_smooth).
## Warning: Removed 20 rows containing missing values (geom_point).

Relationship among testosterone, cortisol, and dyadic sexual desire in women

impact_data_nodrugs_women %>%
  mutate(cortisol_screen_c_q = cut(cortisol_screen_c, breaks = quantile(cortisol_screen_c, probs = seq(0, 1, 1 / 5), na.rm = TRUE), include.lowest = TRUE)) %>% 
  filter(!is.na(cortisol_screen_c_q)) %>% 
  ggplot(mapping = aes(x = testosterone_screen_c, y = SDI_dyad)) +
  geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"), se = FALSE, color = "#CC79A7") +
  geom_point(alpha = 0.25, size = 3, color = "#CC79A7") +
  facet_wrap(facets = ~ cortisol_screen_c_q) +
  scale_y_continuous(breaks = seq(0, 70, 10)) +
  labs(x = "Testosterone (mean-centered)", y = NULL, title = "Sum Dyadic Sexual Desire (6-62 points) by mean-centered cortisol quintile in women") +
  theme(legend.position = "top")
## Warning: Removed 22 rows containing non-finite values (stat_smooth).
## Warning: Removed 22 rows containing missing values (geom_point).

Relationship among testosterone, cortisol, and solitary sexual desire in women

impact_data_nodrugs_women %>%
  mutate(cortisol_screen_c_q = cut(cortisol_screen_c, breaks = quantile(cortisol_screen_c, probs = seq(0, 1, 1 / 5), na.rm = TRUE), include.lowest = TRUE)) %>% 
  filter(!is.na(cortisol_screen_c_q)) %>% 
  ggplot(mapping = aes(x = testosterone_screen_c, y = SDI_sol)) +
  geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"), se = FALSE, color = "#CC79A7") +
  geom_point(alpha = 0.25, size = 3, color = "#CC79A7") +
  facet_wrap(facets = ~ cortisol_screen_c_q) +
  scale_y_continuous(breaks = seq(0, 25, 5)) +
  labs(x = "Testosterone", y = NULL, title = "Sum Solitary Sexual Desire (2-23 points) by mean-centered cortisol quintile in women") +
  theme(legend.position = "top")
## Warning: Removed 17 rows containing non-finite values (stat_smooth).
## Warning: Removed 17 rows containing missing values (geom_point).

Relationship among testosterone, perceived stress, and dyadic sexual desire in men

impact_data_nodrugs_men %>%
  mutate(pss_c_q = cut(pss_c, breaks = quantile(pss_c, probs = seq(0, 1, 1 / 5), na.rm = TRUE), include.lowest = TRUE)) %>% 
  filter(!is.na(pss_c_q)) %>% 
  ggplot(mapping = aes(x = testosterone_screen_c, y = SDI_dyad)) +
  geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"), se = FALSE, color = "#0072B2") +
  geom_point(alpha = 0.25, size = 3, color = "#0072B2") +
  facet_wrap(facets = ~ pss_c_q) +
  scale_y_continuous(breaks = seq(0, 70, 10)) +
  labs(x = "Testosterone (mean-centered)", y = NULL, title = "Sum Dyadic Sexual Desire (6-62 points) by mean-centered perceived stress quintile in men") +
  theme(legend.position = "top")
## Warning: Removed 27 rows containing non-finite values (stat_smooth).
## Warning: Removed 27 rows containing missing values (geom_point).

Relationship among testosterone, perceived stress, and solitary sexual desire in men

impact_data_nodrugs_men %>%
  mutate(pss_c_q = cut(pss_c, breaks = quantile(pss_c, probs = seq(0, 1, 1 / 5), na.rm = TRUE), include.lowest = TRUE)) %>% 
  filter(!is.na(pss_c_q)) %>% 
  ggplot(mapping = aes(x = testosterone_screen_c, y = SDI_sol)) +
  geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"), se = FALSE, color = "#0072B2") +
  geom_point(alpha = 0.25, size = 3, color = "#0072B2") +
  facet_wrap(facets = ~ pss_c_q) +
  scale_y_continuous(breaks = seq(0, 25, 5)) +
  labs(x = "Testosterone (mean-centered)", y = NULL, title = "Sum Solitary Sexual Desire (2-23 points) by mean-centered perceived stress quintile in men") +
  theme(legend.position = "top")
## Warning: Removed 31 rows containing non-finite values (stat_smooth).
## Warning: Removed 31 rows containing missing values (geom_point).

Relationship among testosterone, perceived stress, and dyadic sexual desire in women

impact_data_nodrugs_women %>%
  mutate(pss_c_q = cut(pss_c, breaks = quantile(pss_c, probs = seq(0, 1, 1 / 5), na.rm = TRUE), include.lowest = TRUE)) %>% 
  filter(!is.na(pss_c_q)) %>% 
  ggplot(mapping = aes(x = testosterone_screen_c, y = SDI_dyad)) +
  geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"), se = FALSE, color = "#CC79A7") +
  geom_point(alpha = 0.25, size = 3, color = "#CC79A7") +
  facet_wrap(facets = ~ pss_c_q) +
  scale_y_continuous(breaks = seq(0, 70, 10)) +
  labs(x = "Testosterone (mean-centered)", y = NULL, title = "Sum Dyadic Sexual Desire (6-62 points) by mean-centered perceived stress quintile in women") +
  theme(legend.position = "top")
## Warning: Removed 36 rows containing non-finite values (stat_smooth).
## Warning: Removed 36 rows containing missing values (geom_point).

Relationship among testosterone, perceived stress, and solitary sexual desire in women

impact_data_nodrugs_women %>%
  mutate(pss_c_q = cut(pss_c, breaks = quantile(pss_c, probs = seq(0, 1, 1 / 5), na.rm = TRUE), include.lowest = TRUE)) %>% 
  filter(!is.na(pss_c_q)) %>% 
  ggplot(mapping = aes(x = testosterone_screen_c, y = SDI_sol)) +
  geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"), se = FALSE, color = "#CC79A7") +
  geom_point(alpha = 0.25, size = 3, color = "#CC79A7") +
  facet_wrap(facets = ~ pss_c_q) +
  scale_y_continuous(breaks = seq(0, 25, 5)) +
  labs(x = "Testosterone", y = NULL, title = "Sum Solitary Sexual Desire (2-23 points) by mean-centered perceived stress quintile in women") +
  theme(legend.position = "top")
## Warning: Removed 31 rows containing non-finite values (stat_smooth).
## Warning: Removed 31 rows containing missing values (geom_point).

Models

Fit model 1

lmer_fit1 <- lmer(SDI_dyad ~ bmi_c + time_since_waking_c + relationship_status_fct + illness_fct + Time_AdjZero * (testosterone_screen_c * cortisol_screen_c + testosterone_screen_c * pss_c) + (1 + Time_AdjZero | ID), data = impact_data_nodrugs_men, REML = TRUE)
## Warning: Some predictor variables are on very different scales: consider
## rescaling

## Warning: Some predictor variables are on very different scales: consider
## rescaling

Diagnostic Plots

# residuals vs. fitted
plot(lmer_fit1, type = c("p", "smooth"))

# random effects
dotplot.ranef.mer(ranef(lmer_fit1, condVar = TRUE))
## $ID

Summary

summary(lmer_fit1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: SDI_dyad ~ bmi_c + time_since_waking_c + relationship_status_fct +  
##     illness_fct + Time_AdjZero * (testosterone_screen_c * cortisol_screen_c +  
##     testosterone_screen_c * pss_c) + (1 + Time_AdjZero | ID)
##    Data: impact_data_nodrugs_men
## 
## REML criterion at convergence: 1536.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.29207 -0.52799  0.01013  0.47190  2.90016 
## 
## Random effects:
##  Groups   Name         Variance Std.Dev. Corr 
##  ID       (Intercept)  64.103   8.006         
##           Time_AdjZero  1.003   1.002    -0.31
##  Residual              25.458   5.046         
## Number of obs: 222, groups:  ID, 62
## 
## Fixed effects:
##                                                        Estimate Std. Error
## (Intercept)                                           3.342e+01  3.168e+00
## bmi_c                                                 1.445e-01  2.752e-01
## time_since_waking_c                                  -4.703e-02  1.178e-01
## relationship_status_fct.L                            -4.670e+00  1.187e+00
## relationship_status_fct.Q                            -5.268e+00  1.310e+00
## illness_fct.L                                        -1.079e+01  6.105e+00
## illness_fct.Q                                        -5.470e+00  3.594e+00
## Time_AdjZero                                         -1.598e-01  2.448e-01
## testosterone_screen_c                                -5.509e-02  3.417e-02
## cortisol_screen_c                                     1.113e+00  8.467e-01
## pss_c                                                 1.306e-01  1.109e-01
## testosterone_screen_c:cortisol_screen_c               1.639e-03  1.948e-02
## testosterone_screen_c:pss_c                           4.235e-03  3.645e-03
## Time_AdjZero:testosterone_screen_c                    1.969e-02  7.626e-03
## Time_AdjZero:cortisol_screen_c                       -2.260e-01  1.986e-01
## Time_AdjZero:pss_c                                    2.296e-02  3.068e-02
## Time_AdjZero:testosterone_screen_c:cortisol_screen_c -5.342e-04  3.776e-03
## Time_AdjZero:testosterone_screen_c:pss_c             -3.997e-04  7.930e-04
##                                                              df t value
## (Intercept)                                           9.437e+01  10.549
## bmi_c                                                 5.947e+01   0.525
## time_since_waking_c                                   1.464e+02  -0.399
## relationship_status_fct.L                             1.802e+02  -3.934
## relationship_status_fct.Q                             1.940e+02  -4.022
## illness_fct.L                                         7.226e+01  -1.767
## illness_fct.Q                                         8.020e+01  -1.522
## Time_AdjZero                                          2.053e+01  -0.653
## testosterone_screen_c                                 1.815e+02  -1.612
## cortisol_screen_c                                     1.403e+02   1.315
## pss_c                                                 1.571e+02   1.178
## testosterone_screen_c:cortisol_screen_c               1.748e+02   0.084
## testosterone_screen_c:pss_c                           1.754e+02   1.162
## Time_AdjZero:testosterone_screen_c                    1.283e+02   2.583
## Time_AdjZero:cortisol_screen_c                        1.125e+02  -1.138
## Time_AdjZero:pss_c                                    1.198e+02   0.748
## Time_AdjZero:testosterone_screen_c:cortisol_screen_c  1.225e+02  -0.141
## Time_AdjZero:testosterone_screen_c:pss_c              1.444e+02  -0.504
##                                                      Pr(>|t|)    
## (Intercept)                                           < 2e-16 ***
## bmi_c                                                0.601369    
## time_since_waking_c                                  0.690274    
## relationship_status_fct.L                            0.000119 ***
## relationship_status_fct.Q                            8.27e-05 ***
## illness_fct.L                                        0.081491 .  
## illness_fct.Q                                        0.131902    
## Time_AdjZero                                         0.521315    
## testosterone_screen_c                                0.108676    
## cortisol_screen_c                                    0.190726    
## pss_c                                                0.240658    
## testosterone_screen_c:cortisol_screen_c              0.933045    
## testosterone_screen_c:pss_c                          0.246809    
## Time_AdjZero:testosterone_screen_c                   0.010927 *  
## Time_AdjZero:cortisol_screen_c                       0.257654    
## Time_AdjZero:pss_c                                   0.455768    
## Time_AdjZero:testosterone_screen_c:cortisol_screen_c 0.887734    
## Time_AdjZero:testosterone_screen_c:pss_c             0.615001    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling

Variance Partition Coefficients

lmer_vpc(lmer_fit1)
## # A tibble: 4 x 7
##   grp      var1         var2          vcov  sdcor     vpc vpc_pct
##   <chr>    <chr>        <chr>        <dbl>  <dbl>   <dbl> <chr>  
## 1 ID       (Intercept)  <NA>         64.1   8.01   0.728  72.76% 
## 2 ID       Time_AdjZero <NA>          1.00  1.00   0.0114 1.14%  
## 3 ID       (Intercept)  Time_AdjZero -2.46 -0.307 -0.0280 -2.8%  
## 4 Residual <NA>         <NA>         25.5   5.05   0.289  28.9%

Fit model 2

lmer_fit2 <- lmer(SDI_sol ~ bmi_c + time_since_waking_c + relationship_status_fct + illness_fct + Time_AdjZero * (testosterone_screen_c * cortisol_screen_c + testosterone_screen_c * pss_c) + (1 + Time_AdjZero | ID), data = impact_data_nodrugs_men, REML = TRUE)
## Warning: Some predictor variables are on very different scales: consider
## rescaling

## Warning: Some predictor variables are on very different scales: consider
## rescaling

Diagnostic Plots

# residuals vs. fitted
plot(lmer_fit2, type = c("p", "smooth"))

# random effects
dotplot.ranef.mer(ranef(lmer_fit2, condVar = TRUE))
## $ID

Summary

summary(lmer_fit2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: SDI_sol ~ bmi_c + time_since_waking_c + relationship_status_fct +  
##     illness_fct + Time_AdjZero * (testosterone_screen_c * cortisol_screen_c +  
##     testosterone_screen_c * pss_c) + (1 + Time_AdjZero | ID)
##    Data: impact_data_nodrugs_men
## 
## REML criterion at convergence: 1209.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.02331 -0.46265 -0.06901  0.38740  3.13362 
## 
## Random effects:
##  Groups   Name         Variance Std.Dev. Corr 
##  ID       (Intercept)  18.398   4.2893        
##           Time_AdjZero  0.128   0.3578   -0.41
##  Residual               5.323   2.3071        
## Number of obs: 219, groups:  ID, 62
## 
## Fixed effects:
##                                                        Estimate Std. Error
## (Intercept)                                           7.495e+00  1.378e+00
## bmi_c                                                 2.114e-01  1.411e-01
## time_since_waking_c                                   5.665e-02  5.361e-02
## relationship_status_fct.L                            -4.090e-01  5.508e-01
## relationship_status_fct.Q                            -7.621e-01  6.280e-01
## illness_fct.L                                        -2.830e+00  2.557e+00
## illness_fct.Q                                        -2.565e+00  1.516e+00
## Time_AdjZero                                          1.441e-01  9.663e-02
## testosterone_screen_c                                 1.043e-02  1.586e-02
## cortisol_screen_c                                     7.885e-01  3.873e-01
## pss_c                                                 5.504e-02  5.149e-02
## testosterone_screen_c:cortisol_screen_c               2.787e-03  8.922e-03
## testosterone_screen_c:pss_c                          -2.536e-03  1.742e-03
## Time_AdjZero:testosterone_screen_c                   -2.429e-05  3.347e-03
## Time_AdjZero:cortisol_screen_c                       -1.362e-01  8.792e-02
## Time_AdjZero:pss_c                                   -1.650e-02  1.342e-02
## Time_AdjZero:testosterone_screen_c:cortisol_screen_c -1.008e-03  1.657e-03
## Time_AdjZero:testosterone_screen_c:pss_c              2.663e-05  3.527e-04
##                                                              df t value
## (Intercept)                                           8.781e+01   5.440
## bmi_c                                                 5.896e+01   1.498
## time_since_waking_c                                   1.472e+02   1.057
## relationship_status_fct.L                             1.786e+02  -0.743
## relationship_status_fct.Q                             1.893e+02  -1.213
## illness_fct.L                                         5.942e+01  -1.107
## illness_fct.Q                                         6.851e+01  -1.692
## Time_AdjZero                                          2.703e+01   1.491
## testosterone_screen_c                                 1.730e+02   0.657
## cortisol_screen_c                                     1.401e+02   2.036
## pss_c                                                 1.601e+02   1.069
## testosterone_screen_c:cortisol_screen_c               1.642e+02   0.312
## testosterone_screen_c:pss_c                           1.678e+02  -1.455
## Time_AdjZero:testosterone_screen_c                    1.100e+02  -0.007
## Time_AdjZero:cortisol_screen_c                        9.947e+01  -1.549
## Time_AdjZero:pss_c                                    1.066e+02  -1.230
## Time_AdjZero:testosterone_screen_c:cortisol_screen_c  1.030e+02  -0.608
## Time_AdjZero:testosterone_screen_c:pss_c              1.312e+02   0.076
##                                                      Pr(>|t|)    
## (Intercept)                                          4.76e-07 ***
## bmi_c                                                  0.1395    
## time_since_waking_c                                    0.2924    
## relationship_status_fct.L                              0.4587    
## relationship_status_fct.Q                              0.2265    
## illness_fct.L                                          0.2728    
## illness_fct.Q                                          0.0953 .  
## Time_AdjZero                                           0.1475    
## testosterone_screen_c                                  0.5118    
## cortisol_screen_c                                      0.0436 *  
## pss_c                                                  0.2868    
## testosterone_screen_c:cortisol_screen_c                0.7552    
## testosterone_screen_c:pss_c                            0.1474    
## Time_AdjZero:testosterone_screen_c                     0.9942    
## Time_AdjZero:cortisol_screen_c                         0.1244    
## Time_AdjZero:pss_c                                     0.2216    
## Time_AdjZero:testosterone_screen_c:cortisol_screen_c   0.5443    
## Time_AdjZero:testosterone_screen_c:pss_c               0.9399    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling

Variance Partition Coefficients

lmer_vpc(lmer_fit2)
## # A tibble: 4 x 7
##   grp      var1         var2           vcov  sdcor      vpc vpc_pct
##   <chr>    <chr>        <chr>         <dbl>  <dbl>    <dbl> <chr>  
## 1 ID       (Intercept)  <NA>         18.4    4.29   0.792   79.24% 
## 2 ID       Time_AdjZero <NA>          0.128  0.358  0.00551 0.55%  
## 3 ID       (Intercept)  Time_AdjZero -0.630 -0.411 -0.0271  -2.71% 
## 4 Residual <NA>         <NA>          5.32   2.31   0.229   22.92%

Fit model 3

lmer_fit3 <- lmer(SDI_dyad ~ bmi_c + time_since_waking_c + relationship_status_fct + illness_fct + hormonal_contraceptive_fct + Time_AdjZero * (testosterone_screen_c * cortisol_screen_c + testosterone_screen_c * pss_c) + (1 + Time_AdjZero | ID), data = impact_data_nodrugs_women, REML = TRUE)

Diagnostic Plots

# residuals vs. fitted
plot(lmer_fit3, type = c("p", "smooth"))

# random effects
dotplot.ranef.mer(ranef(lmer_fit3, condVar = TRUE))
## $ID

Summary

summary(lmer_fit3)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: SDI_dyad ~ bmi_c + time_since_waking_c + relationship_status_fct +  
##     illness_fct + hormonal_contraceptive_fct + Time_AdjZero *  
##     (testosterone_screen_c * cortisol_screen_c + testosterone_screen_c *  
##         pss_c) + (1 + Time_AdjZero | ID)
##    Data: impact_data_nodrugs_women
## 
## REML criterion at convergence: 1786.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.60643 -0.43524  0.07023  0.50614  2.61812 
## 
## Random effects:
##  Groups   Name         Variance Std.Dev. Corr 
##  ID       (Intercept)  127.5322 11.2930       
##           Time_AdjZero   0.9864  0.9932  -0.47
##  Residual               18.3164  4.2798       
## Number of obs: 267, groups:  ID, 62
## 
## Fixed effects:
##                                                        Estimate Std. Error
## (Intercept)                                           30.488306   1.806817
## bmi_c                                                  0.654744   0.299189
## time_since_waking_c                                    0.064576   0.091562
## relationship_status_fct.L                             -1.349718   0.814299
## relationship_status_fct.Q                             -4.624644   1.149390
## illness_fct.L                                         -3.853232   1.835839
## illness_fct.Q                                         -2.826055   1.290266
## hormonal_contraceptive_fct1                            1.139215   0.961691
## Time_AdjZero                                           0.102374   0.209092
## testosterone_screen_c                                 -0.338364   0.149392
## cortisol_screen_c                                     -0.498657   0.621110
## pss_c                                                 -0.062762   0.118676
## testosterone_screen_c:cortisol_screen_c                0.083385   0.093879
## testosterone_screen_c:pss_c                           -0.024002   0.020258
## Time_AdjZero:testosterone_screen_c                     0.003839   0.030564
## Time_AdjZero:cortisol_screen_c                         0.259325   0.142796
## Time_AdjZero:pss_c                                    -0.024481   0.024435
## Time_AdjZero:testosterone_screen_c:cortisol_screen_c  -0.018275   0.020608
## Time_AdjZero:testosterone_screen_c:pss_c               0.004724   0.003851
##                                                              df t value
## (Intercept)                                          102.443497  16.874
## bmi_c                                                 62.121612   2.188
## time_since_waking_c                                  172.604791   0.705
## relationship_status_fct.L                            184.799372  -1.658
## relationship_status_fct.Q                            244.096781  -4.024
## illness_fct.L                                        166.995058  -2.099
## illness_fct.Q                                        176.686919  -2.190
## hormonal_contraceptive_fct1                          222.509003   1.185
## Time_AdjZero                                          41.617307   0.490
## testosterone_screen_c                                215.440675  -2.265
## cortisol_screen_c                                    207.821220  -0.803
## pss_c                                                211.625613  -0.529
## testosterone_screen_c:cortisol_screen_c              207.982612   0.888
## testosterone_screen_c:pss_c                          200.254853  -1.185
## Time_AdjZero:testosterone_screen_c                   168.614509   0.126
## Time_AdjZero:cortisol_screen_c                       189.541394   1.816
## Time_AdjZero:pss_c                                   146.635226  -1.002
## Time_AdjZero:testosterone_screen_c:cortisol_screen_c 186.792922  -0.887
## Time_AdjZero:testosterone_screen_c:pss_c             155.218831   1.227
##                                                      Pr(>|t|)    
## (Intercept)                                           < 2e-16 ***
## bmi_c                                                  0.0324 *  
## time_since_waking_c                                    0.4816    
## relationship_status_fct.L                              0.0991 .  
## relationship_status_fct.Q                            7.65e-05 ***
## illness_fct.L                                          0.0373 *  
## illness_fct.Q                                          0.0298 *  
## hormonal_contraceptive_fct1                            0.2374    
## Time_AdjZero                                           0.6270    
## testosterone_screen_c                                  0.0245 *  
## cortisol_screen_c                                      0.4230    
## pss_c                                                  0.5975    
## testosterone_screen_c:cortisol_screen_c                0.3754    
## testosterone_screen_c:pss_c                            0.2375    
## Time_AdjZero:testosterone_screen_c                     0.9002    
## Time_AdjZero:cortisol_screen_c                         0.0709 .  
## Time_AdjZero:pss_c                                     0.3180    
## Time_AdjZero:testosterone_screen_c:cortisol_screen_c   0.3763    
## Time_AdjZero:testosterone_screen_c:pss_c               0.2218    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 19 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

Variance Partition Coefficients

lmer_vpc(lmer_fit3)
## # A tibble: 4 x 7
##   grp      var1         var2            vcov  sdcor      vpc vpc_pct
##   <chr>    <chr>        <chr>          <dbl>  <dbl>    <dbl> <chr>  
## 1 ID       (Intercept)  <NA>         128.    11.3    0.901   90.07% 
## 2 ID       Time_AdjZero <NA>           0.986  0.993  0.00697 0.7%   
## 3 ID       (Intercept)  Time_AdjZero  -5.24  -0.467 -0.0370  -3.7%  
## 4 Residual <NA>         <NA>          18.3    4.28   0.129   12.94%

Fit model 4

lmer_fit4 <- lmer(SDI_sol ~ bmi_c + time_since_waking_c + relationship_status_fct + illness_fct + hormonal_contraceptive_fct + Time_AdjZero * (testosterone_screen_c * cortisol_screen_c + testosterone_screen_c * pss_c) + (1 + Time_AdjZero | ID), data = impact_data_nodrugs_women, REML = TRUE)

Diagnostic Plots

# residuals vs. fitted
plot(lmer_fit4, type = c("p", "smooth"))

# random effects
dotplot.ranef.mer(ranef(lmer_fit4, condVar = TRUE))
## $ID

Summary

summary(lmer_fit4)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: SDI_sol ~ bmi_c + time_since_waking_c + relationship_status_fct +  
##     illness_fct + hormonal_contraceptive_fct + Time_AdjZero *  
##     (testosterone_screen_c * cortisol_screen_c + testosterone_screen_c *  
##         pss_c) + (1 + Time_AdjZero | ID)
##    Data: impact_data_nodrugs_women
## 
## REML criterion at convergence: 1243.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.5128 -0.2872 -0.0430  0.3133  3.0994 
## 
## Random effects:
##  Groups   Name         Variance Std.Dev. Corr
##  ID       (Intercept)  21.74688 4.6634       
##           Time_AdjZero  0.01928 0.1389   0.05
##  Residual               1.78059 1.3344       
## Number of obs: 272, groups:  ID, 62
## 
## Fixed effects:
##                                                        Estimate Std. Error
## (Intercept)                                            5.271079   0.682179
## bmi_c                                                  0.616864   0.131623
## time_since_waking_c                                    0.025806   0.027857
## relationship_status_fct.L                              0.438282   0.252444
## relationship_status_fct.Q                             -0.241520   0.342597
## illness_fct.L                                          0.033198   0.562770
## illness_fct.Q                                          0.148842   0.393304
## hormonal_contraceptive_fct1                           -0.455631   0.318846
## Time_AdjZero                                           0.013051   0.046156
## testosterone_screen_c                                 -0.077775   0.042873
## cortisol_screen_c                                      0.112399   0.176193
## pss_c                                                  0.016828   0.033914
## testosterone_screen_c:cortisol_screen_c                0.048831   0.027259
## testosterone_screen_c:pss_c                           -0.014175   0.005756
## Time_AdjZero:testosterone_screen_c                     0.021953   0.008258
## Time_AdjZero:cortisol_screen_c                        -0.057791   0.038724
## Time_AdjZero:pss_c                                    -0.006305   0.006422
## Time_AdjZero:testosterone_screen_c:cortisol_screen_c  -0.004872   0.005640
## Time_AdjZero:testosterone_screen_c:pss_c               0.003125   0.001024
##                                                              df t value
## (Intercept)                                           88.498836   7.727
## bmi_c                                                 60.100235   4.687
## time_since_waking_c                                  180.743414   0.926
## relationship_status_fct.L                            185.586781   1.736
## relationship_status_fct.Q                            151.603194  -0.705
## illness_fct.L                                        173.719644   0.059
## illness_fct.Q                                        188.885413   0.378
## hormonal_contraceptive_fct1                          239.213891  -1.429
## Time_AdjZero                                          18.438617   0.283
## testosterone_screen_c                                177.764141  -1.814
## cortisol_screen_c                                    164.053622   0.638
## pss_c                                                171.475846   0.496
## testosterone_screen_c:cortisol_screen_c              182.800334   1.791
## testosterone_screen_c:pss_c                          175.413774  -2.463
## Time_AdjZero:testosterone_screen_c                    88.288384   2.658
## Time_AdjZero:cortisol_screen_c                        75.440968  -1.492
## Time_AdjZero:pss_c                                    57.475074  -0.982
## Time_AdjZero:testosterone_screen_c:cortisol_screen_c  92.153371  -0.864
## Time_AdjZero:testosterone_screen_c:pss_c              78.731656   3.051
##                                                      Pr(>|t|)    
## (Intercept)                                          1.62e-11 ***
## bmi_c                                                1.64e-05 ***
## time_since_waking_c                                   0.35548    
## relationship_status_fct.L                             0.08420 .  
## relationship_status_fct.Q                             0.48191    
## illness_fct.L                                         0.95303    
## illness_fct.Q                                         0.70553    
## hormonal_contraceptive_fct1                           0.15431    
## Time_AdjZero                                          0.78051    
## testosterone_screen_c                                 0.07135 .  
## cortisol_screen_c                                     0.52441    
## pss_c                                                 0.62039    
## testosterone_screen_c:cortisol_screen_c               0.07489 .  
## testosterone_screen_c:pss_c                           0.01476 *  
## Time_AdjZero:testosterone_screen_c                    0.00932 ** 
## Time_AdjZero:cortisol_screen_c                        0.13977    
## Time_AdjZero:pss_c                                    0.33030    
## Time_AdjZero:testosterone_screen_c:cortisol_screen_c  0.38988    
## Time_AdjZero:testosterone_screen_c:pss_c              0.00311 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 19 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

Variance Partition Coefficients

lmer_vpc(lmer_fit4)
## # A tibble: 4 x 7
##   grp      var1         var2            vcov  sdcor      vpc vpc_pct
##   <chr>    <chr>        <chr>          <dbl>  <dbl>    <dbl> <chr>  
## 1 ID       (Intercept)  <NA>         21.7    4.66   0.922    92.23% 
## 2 ID       Time_AdjZero <NA>          0.0193 0.139  0.000818 0.08%  
## 3 ID       (Intercept)  Time_AdjZero  0.0326 0.0503 0.00138  0.14%  
## 4 Residual <NA>         <NA>          1.78   1.33   0.0755   7.55%